	implicit doubleprecision (a-h,o-z)
	dimension r_kr0(7,4,120),r_kb0(7,4,120),r_kr(4,120),r_kb(4,120)
	dimension r0(94),b0(7,4,121),r2(130),b2(7,4,131)
	dimension raxe(120),baxe(4,120)
	character maglevel*2,secter*1,el*2
	common /cr0/r0/cb0/b0/ckr0/r_kr0/ckb0/r_kb0
	pai=3.14159265358979d0
	r_m0=1.6605402d-27
	e=1.60217733d-19
	c=2.997924575d8
      open(11,file='saiif.in',status='old')
	read(11,*)el,a,z,w,h
	close(11)
	rex=3.207d0
	w0=931.478d0
	garma=(w+w0)/w0
	beita=dsqrt(garma**2.0d0-1.0d0)/garma
	v=beita*c
	frev=v/(2.0d0*pai*rex)
	write (6,fmt='(2x,a2,2x,a3,f6.1,2x,a3,f5.1,2x,a3,f7.2,a5,/,
	&2x,a6,f10.6,a3,2x,a2,f5.3/)')el,'A =',a,'Z =',z,'W = ',w,'MeV/A',
     &'Frev =',frev/1.0d6,'Mhz','h=',h
	frad=frev*h
	do i=1,7
	do j=1,4
	if (i.eq.1) maglevel='06'
	if (i.eq.2) maglevel='10'
	if (i.eq.3) maglevel='12'
	if (i.eq.4) maglevel='14'
	if (i.eq.5) maglevel='15'
	if (i.eq.6) maglevel='16'
	if (i.eq.7) maglevel='17'
	if (j.eq.1) secter='1'
	if (j.eq.2) secter='2'
	if (j.eq.3) secter='3'
	if (j.eq.4) secter='4'
	open (unit=11,
	&file='kbkr\k'//maglevel//'s'//secter//'p.dat',
	&status='old')
	read (11,*) p1,p2,p3
	read (11,fmt='(10(f8.5))') (r_kr0(i,j,k),k=1,120)
	read (11,fmt='(10(f8.5))') (r_kb0(i,j,k),k=1,120)  
	close (11)
	open (unit=11,file='magmap\b'//maglevel//'a360p.dat',
	&status='old')
	read (11,*)
10	read (11,fmt='(5(f10.3))') (b0(i,j,k),k=1,95)
	if (b0(i,j,95).ne.90.0d0*(j-1)+45.0d0) goto 10
	close (11)
	end do
	end do
	do k=1,35
	r0(k)=84.0d0+2.0d0*dble(k-1)
	end do
	do k=36,71
	r0(k)=152.0d0+4.0d0*dble(k-35)
	end do
	do k=72,94
	r0(k)=296.0d0+2.0d0*dble(k-71)
	end do
	do k=1,130
	r2(k)=84.0d0+2.0d0*dble(k-1)
	end do
	do i=1,7
	do j=1,4
	do k=1,130
	call larg1dr4bf(i,j,r2(k),b_res)
	b2(i,j,k)=b_res
	end do
	b2(i,j,131)=b0(i,j,95)
	end do
	end do
	do i=1,7
	do j=1,4
	do k=1,120
	b0(i,j,k)=b2(i,j,k+(1000-840)/20)
	end do
	b0(i,j,121)=b2(i,j,131)
	end do
	end do
	do k=1,120
	raxe(k)=1000.0d-3+20.0d-3*dble(k-1)
	end do
	do j=1,4
	do k=1,120
	n=0
	i=2
	r_kr1=r_kr0(i,j,k)
	r_kb1=r_kb0(i,j,k)
	baxe(j,k)=b0(i,j,k)
20	n=n+1
	r_kr(j,k)=r_kr1
	r_kb(j,k)=r_kb1
	baxec=baxe(j,k)
	v=2.0d0*pai*frev*raxe(k)/r_kr(j,k)
	garma=1.0d0/dsqrt(1.0d0-(v/c)**2)
	baxe(j,k)=(a/(z*e))*r_m0*garma*2.0d0*pai*frev*r_kb(j,k)
	baxe(j,k)=1.0d4*baxe(j,k)
	write (6,*) j,k,n,baxe(j,k)
	call larg1dr4bl(j,k,baxe(j,k),r_kr_res,r_kb_res)
	r_kr1=r_kr_res
	r_kb1=r_kb_res
	if (n.eq.1) goto 20
	if (dabs(baxe(j,k)-baxec) .ge. 1.0d-7) goto 20
	end do
	end do
	open (unit=12,
	&file='saiifout.dat',status='unknown')
	write (12,fmt=
	&'(2x,a2,2x,a3,f6.1,2x,a3,f5.1,2x,a3,f7.2,a5,/,
	&2x,a4,f10.6,a3,2x,a3,f4.1,/)') 
	&el,'A =',a,'Z =',z,'W =',w,'MeV/A',
	&'RF =',frad/1.0d6,'Mhz','H =',h
	do j=1,4
	write (12,fmt='(1x,i1)') j
	write (12,fmt='(5(f10.3))') (baxe(j,k),k=1,120)
	end do
	close (12)
	open (unit=12,file='saiifrev00.dat',status='unknown')
	write (12,fmt=
	&'(2x,a2,2x,a3,f6.1,2x,a3,f5.1,2x,a3,f7.2,a5,/,
	&2x,a4,f10.6,a3,2x,a3,f4.1,/)') 
	&el,'A =',a,'Z =',z,'W =',w,'MeV/A',
	&'RF =',frad/1.0d6,'Mhz','H =',h
	do j=1,4
	write (12,fmt='(1x,i1)') j
	write (12,fmt='(5(f10.3))') (baxe(j,k),k=1,120)
	end do
	close (12)
	n_it=0
	open (unit=12,
	&file='n_it.dat',status='unknown')
	write (12,fmt='(2x,i2)') n_it
	close (12)
	stop
	end
		
	subroutine larg1dr4bf(l,n,x,b)
	implicit real(8) (a-h,o-z)
	dimension rb0(94),bf0(7,4,121)
	common /cr0/rb0/cb0/bf0
	dimension x0(94),xi(4)
	ix_max=94
	do i=1,ix_max
	x0(i)=rb0(i)
	end do
	do i=1,ix_max-1
	if (x0(i).le.x.and.x.lt.x0(i+1)) ix=i-2
	end do
	if (x.le.x0(1+1)) ix=0
	if (x.ge.x0(ix_max-1)) ix=ix_max-4
	x1=x0(ix+1)
	x2=x0(ix+2)
	x3=x0(ix+3)
	x4=x0(ix+4)
	xi(1)=(x-x2)*(x-x3)*(x-x4)/((x1-x2)*(x1-x3)*(x1-x4))
	xi(2)=(x-x1)*(x-x3)*(x-x4)/((x2-x1)*(x2-x3)*(x2-x4))
	xi(3)=(x-x1)*(x-x2)*(x-x4)/((x3-x1)*(x3-x2)*(x3-x4))
	xi(4)=(x-x1)*(x-x2)*(x-x3)/((x4-x1)*(x4-x2)*(x4-x3))
	b=0.0d0
	do i=1,4
	b=b+xi(i)*bf0(l,n,ix+i)
	end do
	return
	end
	
	subroutine larg1dr4bl(j,k,x,y1,y2)
	implicit real(8) (a-h,o-z)
	dimension bf0(7,4,121),r_kr0(7,4,120),r_kb0(7,4,120)
	common /cb0/bf0/ckr0/r_kr0/ckb0/r_kb0
	dimension x0(7),xi(4)
	ix_max=7
	do i=1,ix_max
	x0(i)=bf0(i,j,k)
	end do
	do i=1,ix_max-1
	if (x0(i).le.x.and.x.lt.x0(i+1)) ix=i-2
	end do
	if (x.le.x0(1+1)) ix=0
	if (x.ge.x0(ix_max-1)) ix=ix_max-4
	x1=x0(ix+1)
	x2=x0(ix+2)
	x3=x0(ix+3)
	x4=x0(ix+4)
	xi(1)=(x-x2)*(x-x3)*(x-x4)/((x1-x2)*(x1-x3)*(x1-x4))
	xi(2)=(x-x1)*(x-x3)*(x-x4)/((x2-x1)*(x2-x3)*(x2-x4))
	xi(3)=(x-x1)*(x-x2)*(x-x4)/((x3-x1)*(x3-x2)*(x3-x4))
	xi(4)=(x-x1)*(x-x2)*(x-x3)/((x4-x1)*(x4-x2)*(x4-x3))
	y1=0.0d0
	y2=0.0d0
	do i=1,4
	y1=y1+xi(i)*r_kr0(ix+i,j,k)
	y2=y2+xi(i)*r_kb0(ix+i,j,k)
	end do
	return
	end
	